Zero density limit extrapolation of the superfluid transition temperature in a unitary atomic Fermi 
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The superfluid transition temperature T c of a unitary Fermi gas on a three-dimensional isotropic lattice with 
an attractive on-site interaction is investigated as a function of density n, from half filling down to 5.0 x 10 -7 
per unit cell, using a pairing fluctuation theory. We show that except at very low densities (n 1 ^ 3 < 0.2), where 
T c /Ep is linear in n 1//3 , T c /Ep exhibits significant higher order nonlinear dependence on n 1 / 3 . Therefore, 
linear extrapolation using results at intermediate densities such as in typical quantum Monte Carlo simulations 
leads to a significant underestimate of the zero density limit of T c /Ep. Our result, T c /Ep — 0.256, at n — 
is subject to reduction from particle-hole fluctuations and incoherent single particle self energy corrections. 

PACS numbers: 03.75.Ss, 03.75.Nt, 74.20.-z, 74.25.Dw 
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Experimental realization of superfluidity in cold atomic 
Fermi gases has given the BCS-Bose-Einstein condensa- 
tion (BEC) crossover study a strong boost over the past 
decade. More importantly, main interests have been paid to 
the strongly interacting regime, where the s-wave scattering 
length a is large. In particular, the unitary limit, where the 
scattering length diverges, has become a test point for theo- 
ries. As a consequence, the superfluid transition temperature 
T c in a unitary Fermi gas has been under intensive investiga- 
tion in recent years. 

Apart from calculating T c directly in the 3D continuum 
with various approximations fll 4l3l . one important method is 
to calculate T c on a lattice and then extrapolate to zero density. 
It has been argued that the zero density limit is identical to the 
continuum case. Indeed, this is the approach used by quantum 
Monte Carlo (QMC) simulations. For this approach to work, 
two conditions have to be met. First, the result obtained from 
the simulation at a given density has to be accurate; this re- 
quires that both the lattice size and the particle number have to 
be large enough. Second, the densities n at which the simula- 
tions are performed have to be in the asymptotic linear regime 
of T c as a function of n 1 ^. 

It is extremely important to investigate this issue, because 
the results of QMC have often been taken with high credibility 
in the cold atom community, despite the large discrepancies 
between the results from different groups (as well as within 
the same group sometimes), and the small total fermion num- 
ber and lattice size used. For example, using QMC, Troyer 
and coworkers 1171 [lal reported T c /Ep — 0.152, whereas 
Bulgac et al. ifl 9L l2C4l reported T C /E F = 0.23 and 0.15 in 
different papers. Using the method of Ref. 17, Goulko and 
Wingate UM found T C /E F = 0.171. Another recent result 
[22] from QMC gave T C /E F = 0.245. Although these dif- 
ferent results do not seem to be converging, the above result 
of Troyer and coworkers lfl7i [l8ll has been widely cited and 
compared with recently. The simulations in Ref.[l7| were done 
for lattice fermions at finite densities and then extrapolated to 
zero density. It is the purpose of the present paper to investi- 
gate how low in density one needs to go so that the simulations 
are in the asymptotic linear regime to ensure the accurateness 



of the zero density limit extrapolation. 

In this paper, we will study the finite density effect on the 
zero density limit extrapolation by calculating T c on a 3D 
isotropic lattice with an attractive on-site interaction, U, using 
a pairing fluctuation theory. This theory has been able to gen- 
erate theoretical results in good agreement with experiment 
lfl2[|23ll . To show how the lattice effect evolves with fermion 
density, we drop the complication of the particle-hole channel. 
Our result reveals that, as the density approaches zero, T c /Ep 
does reach the 3D continuum value. However, linear extrap- 
olation using data points calculated at intermediate densities, 
such as those in Ref. [17], will lead to a significant underesti- 
mate of T c for the continuum limit. When particle-hole chan- 
nel contributions are properly included [24], we expect that 
the zero density limit will yield T c /Ep = 0.217, as directly 
calculated in the continuum. 

Details of the pairing fluctuation theory can be found in 
Ref. [Io| both in the continuum and on a lattice (see Ref. 13 
for the treatment of the particle-hole channel effect). On a lat- 
tice, the fermion dispersion is give by £k = 2i (3 — cosk x — 
cos k y — cos k z )—[i = £k — Mi where t is the hopping integral, 
ek is the kinetic energy, and we have set the lattice constant ao 
to unity. We define Fermi energy for a given density n by the 
chemical potential for a non-interacting Fermi gas at zero T. 
In addition, a contact potential in the continuum now becomes 
an on-site attractive interaction U. Namely, we are now solv- 
ing a negative U Hubbard model. The Lippmann-Schwinger 
relation reads m/AitaTv 2 = l/U + X)k(l/ e k)- Therefore, the 
critical coupling strength is given by U c — — l/Sk^AO = 
— 7.91355t. Here m — t/2 is the effective fermion mass in 
the dilute limit. In what follows, we shall set ks = h= 1. 

To recapitulate our theory, the fermion self energy comes 
from two contributions, associated with the superfluid con- 
densate and finite momentum pairs, respectively, given 



by E(if) 



•Esc(K) + % P9 (K), where X SC {K) 



-A 2 SC G (-K) and E pg (K) = £ Q t pg (Q)G (Q - K), with 
A sc being the superfluid order parameter. Y. SC (K) vanishes 
at and above T c . The finite momentum T-matrix t pg (Q) = 
U/[l + Ux(Q)] derives from summation of ladder diagrams in 
the particle-particle channel, with pair momentum Q, where 
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the pair susceptibility \{Q) = J2k G ( k ) G o{Q - K ) in- 
volves the feedback of the self energy via the full Green's 
function G(K). As usual, we use a four vector notation, 
K = (iw,,k), Q = (ifin.q), Ex = r E/ Ek> and 
Eq = ^En Eq' wnere w z (^n) are the odd (even) Mat- 
subara frequencies. 

By the Thouless criterion, the T c equation, given by 1 + 
f7%(0) = 0, now contains the self energy feedback. This is 
a major difference between our pairing fluctuation theory and 
those based on Nozieres and Schmitt-Rink (NSR) JH or sad- 
dle point approximations |2|] . 

After analytical continuation ifl n — > £1 + i0 + , one can 
Taylor expand the (inverse) T-matrix as <~ 1 (f2, q) ss Z(Vl — 
+ f^pair + *rq)i ar, d thus extract the pair dispersion f2 q = 
2B(3 — cos q x — cos q y — cos q z ). Here the imaginary part T q 
can be neglected when pairs become (meta)stable iflOfl . 

At and below T c , fi pa ir = and H pg (K) can be approxi- 



mated as Y, pg (K) 
with the pseudogap parameter A pg defined as 
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where b(x) is the Bose distribution function. 

Neglecting the incoherent term SH in S pg , we arrive at the 
total self energy T,(K) in the BCS form: 

E(tf) » -A 2 G (-K), (2) 

where the total gap A is determined via A 2 = A 2 C + A 2 g . 
Therefore, the Green's function G{K), the quasiparticle dis- 
persion E k = y£k + anc l tr,e 8 a P *° r ^c) equation all 
follow the BCS form, except that the total gap A now con- 
Joins Z?ofn contributions from the order parameter A sc and 
the pseudogap A pg . Thus the gap equation is given by 



i+uj: 



2K. 



0. 



(3) 



where /(x) is the Fermi distribution function. In addition, the 
number equation, n = 2 Ex G{K), is given by, 



E 



(4) 



Equations ||3), (O, and (fl]i form a closed set. For given 
interaction U, they can be used to solve self consistently for 
T c as well as A and fi at T c . 

In Fig.[T|we plot T c as a function of pairing strength —U /6t 
for various densities from high to low. Here 6t is the half band 
width. For n = 0.7, the maximum T c occurs on the BEC side 
of unitarity. Then it moves to the BCS side as n decreases. 
As n further decreases, the maximum moves slowly back to 
the unitary point. This should be contrasted with the 3D con- 
tinuum case, for which the maximum occurs slightly on the 
BEC side. The fact that the maximum occurs on the BCS side 
manifests strong lattice effect at these intermediate densities; 
it is the lattice effect that causes difficulty for pair hopping and 
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Figure 1 . (Color online) Behavior of T c /6t as a function of the attrac- 
tive on-site interaction —U/6t on a 3D isotropic lattice for various 
density from high to low, as labeled. The unitary limit corresponds 
to —U/6t = 1.31893, as indicated by the vertical dotted line. 



thus suppresses T c . Even at density as low as n — 0.005, the 
maximum is still slightly on the BCS side. 

In order to compare with the continuum T c curves (see 
Fig. 10 in Ref. [l2j for example) more easily, we normalize 
the T c curves by corresponding Fermi energy Ep, as shown 
in Fig. [2] For clarity, we have dropped the curves for the two 
high densities, n = 0.7 and 0.5. The lattice effect has made 
the peak around unitarity much more pronounced, and nec- 
essarily present in all different theoretical treatments of finite 
temperature BCS-BEC crossover lf25ll . As n decreases, this 
peak becomes narrower and moves closer to unitarity. Beyond 
the unitary limit, the curve for n — 0.001, as a low density ex- 
ample, exhibits a rapid falloff with pairing strength, and then 
decreases following the functional form T c oc —t 2 /U. This 
is due to the virtue ionization during pair hopping in the BEC 
regime. At unita rity , a significant fraction of fermions form 
metastable pairs llOlfTlII already at T c , and thus they also see 
the lattice effect during pair hopping through virtue ionization. 
This suggests that the lattice effect will never go away in the 
unitary limit no matter how low the density may be. Figure [2] 
also reveals that, as n approaches zero, the maximum T c / Ep 
as well as T c /Ep at unitarity gradually increase. 

Finally, presented in the main figure of Fig.[3]is T c /Ep as 
a function of (cubic root of) density n in the unitary limit, 
down to n = 5.0 x 10~ 7 , since the lattice effect is expected 
to vary as n 1 ^ 3 to the leading order, namely, T c (n)/ Ep(n) — 
T c (0)/E F (0) -aaon 1 / 3 + o(a 2 /n 2 / 3 ), a is a proportionality 
coefficient. Note that aon 1 / 3 represents the ratio between the 
lattice period and the mean interparticle distance. At the same 
time, T c /Ep and T c /6t are plotted as a function of n in the 
upper inset. It shows T c /Ep increases rapidly near the very 
end of n = 0. The behavior of Ep/Qt is shown in the lower 
inset, in a log-log plot. In units of 6t, both Ep and T c vanish at 
n = and reach a maximum at half filling I2al . In particular, 
Ep = 6t at half filling, as expected. 

Our result reveals that as n decreases from half filling, 
T c /Ep decreases and reaches a minimum of 0.172 around 
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Figure 2. (Color online) T c /Ef as a function of — U/6t on a 3D 
isotropic lattice for various density from n — 0.2 to 0.005. The 
unitary limit corresponds to — U /6t — 1.31893, as indicated by the 
vertical dotted line. 



n = 0.28, and then starts to recover slowly. It does not ac- 
celerate until the very end of n — 0. The main plot suggests 
that T c /Ep eventually does recover its continuum counter- 
part value, 0.256, but the curve exhibits a good linearity only 
for n 1 / 3 < 0.2, i.e. n < 0.008. At n = 5.0 x 10~ 7 , we 
finds T C /E F = 0.254, close to 0.256. Using the data below 
n 1 ' 3 = 0.2, our extrapolation (the green dotted line) leads to 
T C /E F = 0.2557 « 0.256 for the continuum limit. Note that 
the data points for n 1 / 3 > 0.3 shows a rather obvious devia- 
tion from the lower n extrapolation line. This implies that the 
range of density for extrapolation used in Ref. [l7| is still far 
from the asymptotic linear regime. In fact, n 1 / 3 > 0.3 cannot 
be regarded as < 1. Indeed, the recent result of Goulko and 
Wingate 112 111 seems to confirm this point. They pushed their 
simulations down to rt 1 / 3 ps 0.23 (albeit with a big error bar), 
and obtained T c /Ep = 0.173 for the zero density limit us- 
ing a linear extrapolation. One can also see from their Fig. 7 
that, without this lower density data point, they would have 
obtained a lower value for T c /Ep. In addition, their quadratic 
fit would yield T c /Ep js 0.19. Finally, we note that a closer 
look of Fig. 3 of Ref. 17 suggests that the lowest density point 
(also with a big error bar) actually already shows that their 
curve starts to bend upward, away from the straight extrapo- 
lation line. Although not conclusive, this observation agrees 
with the T c /Ep curve in Fig. [3] 

Despite the big difference between our theory and the QMC 
approach, it is reasonable to expect that the lattice effect has a 
rather similar effect on T c /Ep. Therefore, we believe that in 
order to obtain an accurate value of T c /Ep in the zero density 
limit using a linear extrapolation, one needs to perform QMC 
down to n 1 / 3 ~ 0.1 (i.e. n ~ 1.0 x 10~ 3 ) or lower. 

A later paper by Troyer et al. Il8ll claimed that they con- 
firmed their lattice fermion result by working in the continuum 
limit. However, it is likely that the lattice effect was actually 
introduced back through their Eqs. (3) and (4) and the peri- 
odic boundary condition 12711 . Indeed, this has been confirmed 
by Ref. j28ll . which was partly motivated by the preprint l29ll 
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Figure 3. (Color online) T c /Ef as a function of n 1 ^ 3 on a 3D 
isotropic lattice at unitarity. Shown in the lower left inset is Ep/Qt 
as a function of n, and plotted in the upper inset are T c /Ef (black 
solid curve) and T c /6t (red dashed line) as a function of n. The 
(green dotted) linear extrapolation line obtained from fitting using 
data points below n 1/3 = 0.2 yields T C /E F = 0.2557 at n = 0. 



of the current paper. Furthermore, the thermodynamic limit 
value was obtained from a linear fit of only 3 data pointsl In 
this later paper, they performed simulations at unitarity down 
to n « 0.05 (with l = 1), or n 1/3 = 0.37. Unfortunately, this 
density was still far from low enough to enable an accurate ex- 
trapolation. Above this density, the curve in the main figure 
of Fig. [3] shows significant deviation from linearity, caused by 
contributions of order n 2 / 3 and higher. Our result shows that 
the densities used for the QMC simulations in Ref. 01711 were 
not low enough to ensure a good linear zero density limit ex- 
trapolation of T c /Ep as a function of n 1 / 3 . 

The QMC simulations by Bulgac et al. were also done on a 
lattice, and the number of atoms and the lattice sizes were too 
small (e.g. only 50-55 atoms on an 8 3 lattice in Refs. [l9l and 
[20j equivalent to n = 0.1 ~ 0.11 or n 1/3 = 0.46 ~ 0.48) to 
study the dilute limit as we have done here. In the context of 
dynamical mean field theory, Privitera et al lf30h studied the 
importance of nonuniversal finite-density corrections to the 
unitary limit and found that "densities around n ~ 0.05 — 0.01 
are not representative of the dilute regime". 

It should be noted that our Ep is the actual Fermi energy 
in the 3D lattice, whereas Troyer and coworkers [17] simply 
defined Ep — tkp = t(2>im) 2 / 3 . This definition will become 
the true Fermi energy only in the dilute limit. 

Finally, it is interesting to note that on the lattice, the inter- 
action at unitarity, U c , is density independent so that U c /Ep 
will scale to infinity as n approaches 0. In contrast, in the 
continuum, a contact potential can be regarded as the cut- 
off momentum fco — > oo limit of an s-wave interaction, 
U{k) — U6(k a — k), which has U c = — 2tt 2 /mk a . Appar- 
ently, this U c is scaled down to for a contact potential. This 
dramatic contrast for U c between 3D lattice and 3D contin- 
uum seems to suggest that the 3D continuum cannot be sim- 
ply taken as the zero density limit of a 3D lattice. Indeed, the 
fermions are always subject to the lattice periodicity no matter 



4 



how low the density is. 

Without including the self energy feedback in the T c equa- 
tion HH, the NSR theory HIHEd and the saddle point ap- 
proxmation [2] predicted T c /Ep = 0.22. Other approaches 
reported T C /E F « 0.26 [3], 0.15 H, and 0.16 ffl, the last 
of which exhibits unphysical non-monotonic first-order-like 
behavior in entropy S(T). Floerchinger et al. [7] found 
T c /Ep = 0.264 even after including particle-hole fluctua- 
tions. Within the present theory, we reported T c /Ep = 0.256 
EMI- 

Experimentally, the Duke group 13211 . in collaboration with 
Chen et al., found T c /Ep = 0.27 through a thermodynamic 
measurement in a unitary 6 Li gas. Later, they ll33l IMI ob- 
tained 0.29 and 0.21 by fitting entropy and specific heat data 
with different formulas. The latter value was obtained as- 
suming a specific heat jump at T c , which may not be justi- 
fied in the presence of a strong pseudogap at T c (See, e.g., 
Refs. [TH [3 and [36h - According to our calculations, T c at 
unitarity in the trap is only slightly higher than its homoge- 
neous counterpart, 0.272 versus 0.256. Similar small differ- 
ence in T c between trap and homogeneous cases is expected 
from other theories as well. Therefore, these measurements 



imply that the homogeneous T c j Ep is about 0.25 ~ 0. 19. Re- 
cently, Ku et al. 113711 reported T c /Ep w 0.167 for a homo- 
geneous Fermi gas by identifying the lambda-like transition 
temperature. 

Our result demonstrates that the n — > limit of the lattice 
T c / Ep does approach that calculated directly in the contin- 
uum. We expect this to remain true when the particle-hole 
channel contributions are properly included. In that case, we 
obtain T C J Ep = 0.217 at unitarity [24], consistent with some 
of the above experimental measurements. Finally, we note 
that inclusion of the incoherent self energy ST, in our calcula- 
tions would further reduce the value of T c /Ep |24J], bringing 
it closer to the result of Ku et al. H3]. 
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